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Recent advances in ultra-cold atomic Fermi gases make it possible to achieve a fermionic superfluid 
with multiple spin components. In this context, any mean-field description is expected to fail, owing 
to the presence of tightly bound clusters or molecules that consist of more than two particles. Here 
we present a detailed study of a strongly interacting multi-component Fermi gas in a highly elongated 
or quasi-one-dimensional harmonic trap, which could be readily obtained in experiment. By using 
the exact Bethe ansatz solution and a local density approximation treatment of the harmonic trap, 
we investigate the equation of state of the multi-component Fermi gas in both a homogeneous and 
trapped environment, as well as the density profiles and low-energy collective modes. The binding 
energy of multi-component bound clusters is also given. We show that there is a peak in the 
collective mode frequency at the critical density for a deconfining transition to a many-body state 
that is analogous to the quark color superconductor state expected in neutron stars. 
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I. INTRODUCTION 

Recent experimental progress achieved in trapped 
ultra-cold atomic gases provides a great opportunity for 
exploring strongly interacting many-body physics. Ow- 
ing to the molecular Feshbach resonance |l| , the strength 
of the interactions between atoms in different hyper- 
fine states (or species) can be arbitrarily changed from 
strong to weak coupling in a well-controlled manner. 
Moreover, interactions can be effectively tuned by us- 
ing optical lattices with a varying tunneling barrier 0|. 
Consequently, a number of many-body models in con- 
densed matter physics and nuclear physics can be read- 
ily accessed in the ultra-cold atomic gases. In the past 
few years, considerable interest has been focused on 
trapped two-component Fermi gases close to a broad 
Feshbach resonance. In particular, the crossover from 
Bardeen-Cooper-Schrieffer (BCS) fermionic superfiuidity 
of Cooper pairs to Bose-Einstein condensation (BEC) of 
tightly bounded molecules has been explored in great de- 
tail 0, fl H S 0, i, i, M, [m, E, E M, EE, H M [ii • 

Multi-component Fermi gases with more than two 
species can be easily trapped and manipulated as well. 
For this type of Fermi gas, bound multi-body clusters are 
expected to appear above a critical interaction strength 
[19(1, in addition to the formation of Cooper pairs due 
to two-body "pairing correlations". As is well known, 
a landmark theoretical result in quantum physics is Efi- 
mov's prediction of a universal set of bound trimer states 
appearing for three identical bosons with a resonant two- 
body interaction jl^l- The existence of such an Efimov 
resonance greatly changes the properties of dilute Bose 
observed recently in an ultracold gas of cesium 
atoms [2l[. Similarly, multi-body clusters are of funda- 
mental importance to fermionic superfluidity in multi- 
component Fermi gases. While Cooper pairs are dom- 



inant in the weak coupHng Hmit, an exotic superfluid 
state with bound multi-body clusters should emerge in 
the strong coupling regime. In between, a quantum phase 
transition is then anticipated to take place (as specifled 
for the three component case in Ref. [22]), in contrast 
to the smooth crossover observed in the two-component 
case [s^]. 

This issue is relevant to outstanding problems in nu- 
clear and particle physics. The quark model of nuclear 
matter at low density describes nucleons as three fermion 
clusters: tri-quark bound states. At sufficiently high den- 
sity and pressure, it is conjectured that a phase-transition 
occurs to a deconfined color superfiuid phase of quark 
matter [13, [131 ■ This is beheved to occur in the inte- 
rior of neutron stars and possibly in heavy- ion coUisions. 
While quark matter has many other features, it is inter- 
esting to find that a physically accessible system of in- 
teracting ultra-cold fermions is also expected to display 
these features of multi-body fermion clusters and quan- 
tum phase-transitions. Ultra-cold atoms could provide a 
means to test these theoretical predictions. If confirmed, 
this would lend support to current ideas in particle the- 
ory and astrophysics. 

The description of multi-body bound states is beyond 
the generally adopted mean-field framework, which treats 
competition between two-body correlations among differ- 
ent atom pairs. This, however, is of importance only in 
the weak coupling Hmit. For strong coupling, the use of 
numerical Monte Carlo techniques is hampered by the 
fermionic sign problem f25l|. It is therefore of great im- 
portance to have an analytically soluble model of multi- 
component Fermi gases, and to study both the weakly 
and strongly interacting regimes on an equal footing. 
An exact analysis is provided in this paper for the case 
of multi-component Fermi gases in one dimension (ID), 
where multi-body bound clusters are always present, re- 
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gardless of the interaction strength. Although the one di- 
mensional analysis gives a smooth crossover rather than 
a true phase transition, we believe that considerable in- 
sight may be obtained for a three dimensional gas as well. 

As well as being exactly soluble, the ID problem has 
great experimental relevance. A quantum degenerate 
trapped two-component Fermi gas in quasi ID has been 
demonstrated recently by loading an ultra-cold Fermi gas 
into a two-dimensional optical lattice [IHl of trapping 
'tubes'. In this configuration, the radial motion of atoms 
perpendicular to each tube is frozen to zero-point oscil- 
lations due to tight transverse confinement, while axial 
motion is only weakly confined. One then obtains an ar- 
ray of effective ID systems, each in an axial harmonic 
trap. The manipulation of more than two species in ID 
is within reach of present-day technology, and is likely to 
be achieved soon in experiments. 

A typical example of these developments is lithium gas 
[27| . which has favorable collisional properties among its 
lowest three hyperfine spin states, denoted |2) and 
1 3), respectively. A recent accurate measurement of the 
scattering lengths between these hyperfine states shows 
that the background interactions are anomalously large 
[27| . with background scattering lengths about — ISOOoo, 
where ao {— 0.0529177 nm) is the Bohr radius. There 
are also three broad s-wave Feshbach resonances located 
at the positions B — 834,811, and 690 Gauss for (1,2), 
(2,3), and (3,1) channels, respectively. These peculiar 
collisional properties are useful to cool the gas down to 
the quantum degenerate regime. Ideally, one expects 
experimentally accessible lowest temperatures for this 
three-state mixture to be in the same range as for two- 
component ensembles, i.e., T ~ O.GSTp, where Tp is 
the Fermi temperature of an ideal Fermi gas. Thus, a 
novel fermionic multi-component superfiuid may be an- 
ticipated. For this reason, three-component lithium gas 
has attracted a great deal of theoretical interest, includ- 
ing analysis of mean-field stat es 1281. l29l. ISOl. 31 , 132] as 
well as phase diagrams [13, [sl, [IJlii,!^, 113, H [39] . 

Here, we report on properties of a multi-component 
Fermi gas in ID. Firstly, using the exact Bethe ansatz 
solution [13, El, m, [i^l, we investigate the exact ground 
state of a homogeneous gas with attractive interactions 
at zero temperature. To make contact with experiments, 
we then consider an inhomogeneous Fermi gas under har- 
monic confinement, within the framework of the local 
density approximation (LDA). The equation of state of 
the system in both the uniform and trapped case are 
investigated in detail. Particular attention is drawn to 
the density profiles and low-lying collective modes of the 
trapped cloud, which are readily measurable in experi- 
ment. We show that the gas becomes more attractive 
as the number of species increases, demonstrating the 
strongly interacting nature of multi-body bound clusters. 

A ID multi-component Fermi gas in an optical lattice 
was considered recently by Capponi and co-workers [37| . 
using both the analytic bosonization approach and the 
numerical density matrix renormalization group method. 



This lattice version of the ID Fermi gas resembles the 
system we consider here. In particular, it also features a 
molecular superfiuid phase in the low density Hmit with 
a strongly attractive interaction [37j . 

The paper is organized as follows. In the following 
section, we outline the theoretical model for a ID multi- 
component Fermi gas. Of particular relevance for an ex- 
perimental reaHzation is our calculation of the effective 
ID coupling constant for the three-component lithium 
gases. In Sec. Ill, we present the exact Bethe ansatz 
solution and discuss the equation of state and the sound 
velocity of a uniform system at zero temperature. In Sec. 
IV, using the LDA method we investigate the density pro- 
file and the equation of state in the trapped environment. 
Also, we describe the dynamics of trapped gases in terms 
of ID hydrodynamic equations and develop a novel algo- 
rithm to solve these equations. The behavior of low-lying 
collective modes is then obtained and discussed. We end 
with some concluding remarks in Sec. V. An appendix 
is used to outline the details of the algorithm used in 
solving the ID hydrodynamic equations. 



II. MODELS 

A quasi-lD multi-component Fermi gas in a highly 
elongated trap can be formed using a two-dimensional op- 
tical lattice [2h|. By suitably tuning the lattice depth, the 
anisotropy aspect ratio X — ujz/ojp oi two harmonic fre- 
quencies can become extremely small. This means that 
the Fermi energy associated with the longitudinal motion 
of the atoms is much smaller than the energy level sepa- 
ration along the transverse direction, i.e., Nhuj^ <C hujp 
and ksT <^ fkup, where N is the total number of atoms. 
In this limit, the transverse motion will be essentially 
frozen out, and one ends up with a quasi-one dimensional 
system. 



A. Interaction Hamiltonian 

We study a gas with pseudo-spin S — {k ~ 1) /2, where 
K (> 2) is the number of components. From now we shall 
assume that the fermions in different spin states attract 
each other via the same short-range potential giDS{x). 
Denoting the mass of each fermion as to, with a total 
fermion number N — J2'i=i (where Ni is the num- 
ber of fermions with pseudo-spin projection I) the first 
quantized Hamiltonian for the system is therefore 
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represents the part of Hamiltonian in free space with- 
out the trapping potential m.uP'x^ jl, while w = is an 
oscillation frequency in the axial direction. There is an 
inter-particle attraction between any two fermions with 
different quantum numbers. 

In an elongated trap, the ID effective coupling constant 
g\D is related to the 3D scattering length a-^o- It is shown 
by Bergeman et al. [ii, [i^ that 
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where Up = yjhj (mujp) is the characteristic oscilla- 
tor length in the transverse axis. The constant A = 
-C(l/2)/%/2 ~ 1.0326 is responsible for a confinement 
induced Feshbach resonance [40| , which changes the scat- 
tering properties dramatically when the 3D scattering 
length is comparable to the transverse oscillator length. 
It is convenient to express gm in terms of an effectivelD 
scattering length, gio = —2h?/ {maio), where 
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Note that in this definition of the ID scattering length, 
our sign convention is opposite to the 3D case. 

In the homogeneous case, we measure the interactions 
by a dimensionless coupling constant 7, which is the ratio 
of the interaction energy density eint to the kinetic energy 
density etin In the weak coupling, Cmt ~ gion and 
efcm ~ h^k'^/{2m) ~ h?n?/m, where n is the total linear 
density of the gas. Therefore, one finds 
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Thus, 7 <C 1 corresponds to the weakly interacting limit, 
while the strong coupling regime is realized when 7^1. 



B. Cluster states 

In the case where all the fermions are distinct - which 
is only possible if the number of fermions is less than or 
equal to k - the spatial wave-function can be completely 
symmetric. This allows one to construct eigenstates with 
identical symmetry to the exact solutions already known 
for a one-dimensional Bose gas with attractive interac- 
tions [47]. This exceptionally simple limiting case gives 
useful physical insight into the multi-particle clusters. 
These will be an essential feature in the physical prop- 
erties of more general solutions. Accordingly, we may 
consider as a trial solution the localized quantum soliton 
state with wavefunction: 
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On calculating the effect of the many-body Hamilto- 
nian, we find that: 
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where the energy £"„ is: 
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This symmetric state with an asymptotic exponentially 
decaying wavefunction in each coordinate direction is the 
fundamental bound cluster. The requirement for this to 
be an eigenstate is simply: 
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In terms of this characteristic length-scale of aio — 1/c, 
one finds that the fundamental cluster binding energy 
can be written as: 



h'^Hi{Hi'^ - 1) 



(2.10) 



These bound clusters can have either a fermionic or 
bosonic character, depending on whether k is odd or 
even. The binding energy scales quadratically with the 
Hamiltonian coupling, and cubically with the number of 
bound fermions, k. Clusters are localized relative to the 
center of mass, with a characteristic length scale of am- 
One may reasonably expect some kind of physical tran- 
sition to occur when the linear density exceeds l/anj, 
since at high density the Pauli exclusion principle will 
not allow independent clusters to form. 

Thus, we can expect this type of symmetric bound 
state to predominate at low density, with a transition 
to a radically different behaviour at high density. This 
transition due to the Pauli exclusion principle is a unique 
feature of a ID attractive Fermi gas, and does not occur 
in a ID attractive Bose gas. By contrast, in an attrac- 
tive Bose gas, clusters or quantum solitons can form with 
arbitrary particle number. They have been already ob- 
served with up to for = 10'' for photons in optical 
fibers [4^, and N = 10^ in the case of ultra-cold atoms 
[49| . However, clearly this is not to be expected in the 
case of fermions, where the spin multiplicity limits the 
size of this type of symmetric cluster. 

For K = 3, here is a close analogy between the sym- 
metry properties of these bound clusters and the color 
symmetry properties of quark models in particle physics. 
In the case of quark matter, free nucleons are three-quark 
bound states. However, at high density, it is conjectured 
that there is a quantum phase-transition to a deconfined 
color superconductor state [13, US]- This is expected 
physically at the core of massive neutron stars or in 
heavy- ion collisions. Although we are interested mainly 
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in the one-dimensional case, where mostly a crossover 
would be expected due to dimensionality, we will show 
that a similar type of deconfining "transition" occurs here 
also. 



C. Harmonic Trap 

In the presence of a harmonic trap, we may charac- 
terize the interactions using the dimensionless coupling 
constant at the trap center 70. For an ideal Fermi gas 
with equal spin population in each component, the total 
linear density is 
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in the large- iV Thomas-Fermi (TF) limit, where 
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are the center linear d ensity an d the TF radius respec- 
tively . Here aho = \/h/{muj) is the characteristic os- 
cillator length in the axial direction. We therefore find 
that: 
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To remove the dependence on the number of components 
K, we define a dimensionless parameter 



S = 



N- 
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to describe the interactions. Note that the parameter 
5 depends inversely on the total number of particles. 
Hence, somewhat counter-intuitively, the gas becomes in- 
creasingly non-ideal with a decreasing number of atoms. 
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Figure 1: (Color online) One-dimensional interaction param- 
eters for a three-component lithium gas in a two-dimensional 
optical lattice above the Feshbach resonances. Above a mag- 
netic field B = 1000 Gauss, the difference of interaction pa- 
rameters between different channels becomes very small, i.e., 
[5i2 - S23)/5i2 < 0.06 and (5i2 - (5i3)/5i2 < 0.16. As a result, 
a three-component lithium gas can be well described by our 
exactly soluble model of ID Fermi gases with a symmetric 
interaction. 



Detailed values of the background scattering length ab, 
resonance position Bq, resonance width W , and leading- 
order correction a have been measured precisely by 
Bartenstein et al. [l^l for all three channels. As an exam- 
ple, we take a trapping frequency of w = = 27r x 400 
Hz, which gives rise to aho = 2052 nm. The number of 
atoms in one tube of the lattice can be approximately 
of order N ~ 100. Given these parameters, we use Eqs. 
(|2.4p and l|2.15p to calculate 8. The estimated dimension- 
less coupHng constants are shown in Fig. (U). Here, we 
focus on the BCS side of the resonances, and thus take 
a magnetic field B > 834 Gauss. We find that S 0(1), 
i.e., the gas is in an intermediate interacting regime. The 
difference in the interactions between different channels 
turns out be small above 1000 Gauss. This justifies our 
choice of the same contact interaction potential between 
different hyperfine spin states. 



D. Parameter values 



For experimental relevance, we now estimate inter- 
action parameters for the on-going experiments on ID 
lithium gases. Consider a gas of ^Li atoms in 2D optical 
lattice with a typical lattice spacing periodicity d — 532 
nm. The transverse oscillator length Op is related to d 
via ttp = d/{-iTS^/'^) |50], where s is the ratio of the lattice 
depth to the recoil energy. Taking s — 10, this equation 
yields ftp ~ 95 nm. Empirically, the 3D scattering length 
of ^Li gases at these broad resonances is given by 

asrf = afc[l + W/{B - B^)][\ + a{B - Bo)] . (2.16) 



III. HOMOGENEOUS MULTI-COMPONENT 
FERMI GASES 

We first consider a uniform multi-component Fermi gas 
in one dimension with symmetric inter-component inter- 
actions. In this case, the model is exactly soluble via the 
Bethe ansatz 0, \M, \M, S, H M, M, M, El, [H, [H, [III . 
For simplicity, we assume that each component has the 
same number of particles, i.e., Ni = N/k {I — 1,2, ■ ■ ■ , 
K = 2S + 1). 
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A. Ground State 

In the ground state, the particles partition into groups 
of K fermions. In each group, the fermions all have differ- 
ent quantum numbers, and are bound together to form 
a K-body cluster. Introducing a linear number density, 
n = N/L^ where L is the size of the system, the ground 
state energy i?hom in the thermodynamic limit is given 
by Elli, 



hom 
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1— \ dklnk^- 
2m 



-c2 p(A), (3.1) 



where the coupling c — ,^ , / ^ — ^i^iu hk^^j " 
quasi-momentum distribution of K-body clusters with 



717/2 = l/fliD and p(A) is the 



cut-off rapidity Q. The quasi-momentum distribution is 



cut-on rapidity ihe quasi-momentum 
determined by an integral equation ji^, [Hi 
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n = K dAp{A), 
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which fix the value of the cut-off rapidity. The last term 
in iJhom is simply the contribution from K-body bound 
states and is equal to —{n/K)eKb, with binding energy 
identical to that given in the single-cluster result of Eq 
(l2lnll : 
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We note that the binding energy of multi-component 
clusters increases rapidly with an increasing number of 
species n. In particular, when k > 3 it is larger than 
the pairing energy of k (k — 1) /2 Cooper pairs. In other 
words, if the binding energy was solely due to Cooper 
pairing, one would expect ecp = k{k— 1) 626/2, where 
£26 — fi^/maljj is the two-body binding energy. The in- 
crease above this level is due to K-body correlations: k 
particles interact more strongly with each in a cluster 
than as isolated pairs of particles. 



B. Sound velocity 

Once the ground state energy is obtained, we calculate 
the chemical potential ^hom — dEhom/dN and the corre- 
sponding sound velocity Chom — \/n{dfihom/dn)/'m. The 
sound velocity will be utiHzed later for predicting measur- 
able collective mode frequencies. For numerical purposes, 
it is convenient to rewrite these in a dimensionless form 



that depends on the dimensionless coupling constant 7 
only. 
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These are related by, 



P{l) = 3e(7)-7 , 
c(7)^.(7)-i4^. 
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It is easy to see that for an ideal multi-component Fermi 

gas. 
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Therefore, as units of energy and sound velocity, we de- 
fine a Fermi energy Ef.k = [h^n^ / {2m)]{TT'^ / k^) and a 
Fermi velocity i'f,k = [hn / m]{T: / k) . 



C. Numerical solutions 

The integral equation for the quasi-momentum distri- 
bution has to be solved numerically. This was partly 
carried out by Schlottmann in the k = 4 case, but for re- 
stricted values of the linear density n and coupling c [s^ . 
Here, for completeness, we solve the integral equation for 
a general coupHng constant 7. The asymptotic behavior 
in the weak and strong coupling limits, not explored in 
literature so far, will be discussed in detail. 

To make the equation dimensionless, let us change vari- 
ables as follows [4011: 



K = Qx\ 2c=QX; p{A)^g{x). 



(3.13) 



In terms of the new variables the quasi-momentum distri- 
bution, normalization condition, and ground state energy 
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become, respectively, 
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To obtain better numerical accuracy for the chemical po- 
tential, it is useful to calculate the derivative of e (7) 
directly. With this goal, we introduce = dX/dj and 
gd{x) = dg{x)/d'y, which satisfy the coupled equations 



Ad5(x) -{iXf + {x-xf 
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The derivative of e (7) is then obtained from 
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Numerically, the two set of integral equations, Eqs. 
p.l4p and Eqs. p.lSp . have been solved by decom- 
posing the integrals on a grid with M = 1024 points 
{xi\Xi e [—1, +1]}. For 5(2;), we start from a set of trial 
distributions g^°'{xi), with corresponding parameters of 
A*^°^ Using the standard method for the integrals {40| . 
we obtain a new distribution g{xi), and update A ac- 
cordingly. We iterate this procedure until g{xi) agrees 
with the previous distribution within a certain tolerance, 
and finally calculate the energy function 6(7) using Eq. 
(|3.14p . The integral equation of gd{x) can be solved in 
the same manner, and finally Eq. I|3.16p gives the deriva- 
tive of the energy function de/d'y. We find that these it- 
erative procedures for solving the integral equations are 
very stable. To obtain the sound velocity, the derivative 
of the chemical potential has been calculated accurately 
as a numerical derivative. 

For an illustrative purpose, we show in Fig. ^ the 
quasi-momentum distribution function g{x) as a function 
of the coupling constant (Fig. [2K) or as a function of the 
number of components (Fig. [2b)- As g{x) is an even 
function, we plot only the part with a positive x. For a 
large interaction strength, its approaches k/2tt, while in 
the weak coupHng Hmit, it reduces to l/(27r). Below we 
discuss this limiting behavior in detail. 
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Figure 2: (Color online) Quasi-momentum distributions as a 
function of the dimensionless coupling constant (panel a) or 
as a function of the number of components (panel b). The dis- 
tribution approaches to ii/2n for a large interaction strength. 
While in the weak coupling limit, it goes to l/(27r), but with 
a sharp increase at the boundary of a; = ±1. 



D. Strong coupling limit 

For a strongly interacting or equivalently, a low density 
gas, in which the dimensionless coupHng constant 7^1, 
the value of A in Eq (|3.15p is extremely large. Thus, the 
integral kernel lX/Tr/[(lX)'^ + {x — x')'^] becomes essentially 
a constant, l/(7r£A). In addition, the quasi-momentum 
distribution function g{x) ~ go. Then, the integral equa- 
tion reduces to 
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At the same time, A = 2k75o. Denoting 
(1/k^) J2e=i find that up to the order 1/7'^, 
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Note that the factor decreases rapidly as the number 
of components increases, and goes like Sk, — In k/k^ when 
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Figure 3: Schematic diagram of interacting bound clusters for 
K = 3, in the low density limit. 



K » 1. It is then straightforward to obtain that, 
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E. Weak coupling limit 

The asymptotic behavior in the weak coupUng Hmit is 
more subtle. Numerical calculation in the small 7 limit 
suggests that g{x) l/(27r) and A ~ 7 ^ 0. We then 
expand the quasi-momentum distribution, 



9{x) = — + f{x), 



where / (x) — 0{j) <C 1 satisfies, 

+1 



(3.21) 



K- 1 

27r 



A 



1 



K-l +1 

/J 1 



{exf + {x~ x'f 2n 

— 1 

(tA) + (X — x'j 



and in dimensional form, 



-^hom 

N 



2rn 3k4 



Cho 



K 

hn TT 



2m k4 } 

7 7 



45^ 

7 

165. 

37 
1/2 



1252 



7" 
2052 



7" 



(3.20) 



The first term on the right hand side of the total en- 
ergy and chemical potential is simply the binding energy 
per particle of a ID bound cluster from Eq (|2.10p . while 
the second term arises from interactions between clusters. 
These are strongly repulsive due to the Pauli exclusion 
principle between identical fermions. In the infinite cou- 
pling constant limit, the system behaves like a spinless 
Tonks-Girardeau gas of bound clusters with hard-core re- 
pulsive interactions. This is shown schematically in Fig 

Each cluster has a density n/n and mass Km [HHl, giv- 
ing rise to a chemical potential h^i:'^{n/ / {2Km) = 
[h? /2m]{-K'^ / K^). This is exactly k times the second term 
in the chemical potential ^hom- It is worth emphasiz- 
ing that the compressibility of the system remains pos- 
itive for 7 ^ 00 as indicated by the first term in the 
sound velocity Chom, {hn/m){TT/ k'^). This means that a 
ID multi-component Fermi gas is mechanically stable, 
even in the strongly attractive regime. In contrast, the 
mechanical stability of a strongly interacting 3D multi- 
component gas with K > 3 is not known exactly. It may 
experience collapse for sufficient large number of compo- 
nents, as suggested by Heiselberg [59|. The parity of the 
number of species k may also play an important role on 
infiuencing the stability in this limit. For the stable, ID 
case treated here, this low-density regime is analogous to 
the regime of isolated nucleons in QCD. 



As A ^ 0, to the leading order £X/7r/[{£X)'^ + {x - x')"^] 
6{x-x'). Thus, 



fix) 



27rn 



dx' {iX/ir) 



fr^J^ (exy + ix- x'Y 

(3.23) 

For small A, the integral in / (x) is well approximated by. 



-1-1 



dx'- 



TT 7 {IXf + {X- X'f 



1 - 



2 £X 



TT 1 — X^ 



(3.24) 



We find then to the order of 7, 

K-l X 



fix) 



27r2 l-a;2' 



(3.25) 



which diverges at the boundary x — ±1. It is straight- 
forward to show that, 



1 K — 1 

dxgix) = — ^AlnAH (3.26) 

2 K-l 



+1 



dx {1 - x'^) g (x) = — 



-A + 



(3.27) 



Substituting these results into the energy function 
6(7), it is easy to obtain. 



^(7) = :^-^^A-i^A2ln2A. (3.28) 



3k2 4^2 
Recall that X — k^ J^i 9 ix) ~ kj/tt. The equation of 
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Figure 4: (Color online) Dependence of the uniform ground 
state energy per particle (a), the chemical potential (b), and 
the velocity of sound (c) on the dimensionless coupling con- 
stant 7, at several number of species as indicated. The en- 
ergy and sound velocity are in units of the Fermi energy 
eF,K = [fi^n^/(2m)](7r^/fi:^) and the Fermi velocity vf,k = 
[hn/m]{n / k), respectively. Thin solid lines are the analytic 
results in the two limiting cases for ft = 2, as described in Eq. 
Il3:20ll and Eq. p^Mjl . 



state and the sound velocity are finally given by, 
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(3.29) 



where the first terms on the right hand side are identi- 
cal to an ideal (non-interacting) multi-component Fermi 



gas, as one might expect. Note that the binding energy 
of bound clusters is of second order in 7, and therefore 
is not included in the above expressions. Two-body cor- 
relations are dominant in the weakly interacting limit, 
and give rise to the mean-field Hartree-Fock attractive 
corrections in the second term on the right hand side. 
The non-perturbative terms of order 7^ In^ 7 are beyond 
mean-field theory. This regime is analogous to the color- 
superconductor regime expected in quark matter. 



F. Numerical results 

In Fig. llH), we give the equation of state and the 
sound velocity as a function of the dimensionless cou- 
pling constant 7, obtained by numerically solving the in- 
tegral equations. The ground state energy per particle 
and the chemical potential are measured in units of one- 
third of the Fermi energy and the Fermi energy, respec- 
tively, while the sound velocity is in units of the Fermi 
velocity. We consider a number of components ranging 
up to K = 5 to show the overall trend. 

Starting from the ideal gas results, the thermodynamic 
and dynamic quantities decrease with increasing coupHng 
constant, and finally saturate to the Tonks-Girardeau 
(repulsive) gas limit, as already anticipated. The rate 
of decrease is much faster as the number of components 
K increases. This impHes that the gas becomes more at- 
tractive when the number of particles in the bound clus- 
ter increases. We show also in the figure the asymptotic 
behavior in the two Hmiting cases for k = 2 by thin soHd 
lines. These fit fairly well with the full numerical results, 
apart from a small intermediate interaction region about 
7-1. 

It is worth noting that in the strongly interacting limit, 
the properties of the uniform gas do not exhibit any 
(even-odd) parity dependence on the number of species 
K, because of Tonks-Girardeau fermionization, as men- 
tioned earlier. For a ID multi-component gas in an op- 
tical lattice, however, an entire different picture emerges 
[131 ■ Due to the localization of the lattice in the strong at- 
tractive regime, the system either becomes an ideal Fermi 
gas for odd k, or an ideal Bose gas for even k, exhibiting 
a distinct parity effect. 



IV. MULTI-COMPONENT TRAPPED FERMI 
GAS 

To make quantitative contact with on-going experi- 
ments, it is crucial to take into account the external har- 
monic trapping potential Vtrap{x) = muj'^x^ /2 , which is 
necessary to prevent the fermions from escaping. For a 
large number of fermions, which is likely to be ~ 100 
experimentally, an efficient way to take the trap into ac- 
count is by using the local density approximation (LD A) . 
Together with the exact homogeneous equation of state 
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of a ID multi-component Fermi gas, this gives an asymp- 
totically exact results as long as 1. 

The basic idea of the LDA is that an inhomogeneous 
gas of large size can be treated locally as infinite matter 
with a local chemical potential. We may then partition 
the cloud into many blocks, in each of which the num- 
ber of fermions is much greater than unity. Provided 
the variation of the trap potential across each block is 
negligible compared with the local Fermi energy, any in- 
terface effects may be safely neglected. Thus, each block 
is un-correlated with the others. We note that in ID the 
interface energy scales as compared to the total en- 
ergy, and thus the LDA becomes valid provided iV ^ 1. 

In detail, the LDA amounts to determining the chem- 
ical potential /i from the local equilibrium condition 



where 5 is the interaction parameter for a trapped gas 
defined earher in Eq. I|2.15p . It is clear that the LDA 
equations are controlled by a single parameter (S: (5 <C 1 
corresponds to the weakly coupling limit, while S ^ I 
corresponds to the strongly interacting regime. 

The numerical procedure of solving the LDA equations 
is straightforward. For a given interaction parameter S, 
and initial guess for Jig , we invert the dimensionless local 
equilibrium equations to find j{x) and the linear density 
n{x) — 2/7 (i). The chemical potentials Jig are then 
adjusted to enforce the number conservation, giving a 
better estimate for the next iterative step. The iteration 
is continued until the number normalization condition is 
satisfied within a certain range. 



A. Density profiles and the equation of state 



Mhom [nix)] + -muj'^x'^ = ^g, 
under the normalization restriction. 



N 



n (x) dx, 



(4.1) 



(4.2) 



where n (x) is the total Hnear number density and is 
nonzero inside a Thomas- Fermi radius xtf- We have 
used the subscript "g" to distinguish the global chemi- 
cal potential /ig from the local chemical potential /ihom. 
Rewriting /ihom into the dimensionless form ^[j{x)], 
where 7 (x) — 2/[n (x) aic], we find that 



(k2 _ 1) fi2 ^2^2 



1 



2m t^hix)] + ^muj^x^ ^ Hg. (4.3) 



The first term on the left hand side is simply the bind- 
ing energy, and causes a constant shift to the chemical 
potential. To solve the LDA equations, it is simplest to 
transform into a dimensionless form, by defining 



'■ID 



fi2 

naiD, 



(4.4) 
(4.5) 
(4.6) 



where the binding energy is now absorbed in the re- 
definition of chemical potential. Thus, the local equi- 
librium condition becomes. 



n'^{x) 



H [7 (x)] 



x^ 



Ms: 



(4.7) 



where the dimensionless coupling constant now takes the 
form, 7 (i) — 2/n{x). Accordingly, the normalization 
condition is given by 



+ XTF 



dih(x) = 

at 



1 

J2' 



(4.; 



The asymptotic behavior of density profiles can be de- 
termined analytically in the strong and weak coupHng 
limits. In the strong interaction regime of ^ ^ 1, 
^((7) = (7rVK^)[l + 16S'«/(37)], and the local equilib- 
rium condition is given by. 



n (x) TT 



1 



37 (i) 



Aig- 



(4.9) 



In the infinite coupHng limit of Tonks- Girardeau gas, 
where 7 (i) 00 and therefore the second term in /i (7) 
vanishes, the density profile takes the form. 



flTG (x) 



V2h. 



2r2 



1 
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n6 \ 2 
and the global chemical potential is. 
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(4.10) 



(4.11) 



The inclusion of the next order of I/7 in /i (7) leads to 
a density variation Ah{x), as well as a change in the 
chemical potential AJlg. Linearizing the local equilibrium 
condition, we find that. 



An (x) 



Aflg 



8Sk r_ ._.i2 

— [UTG [X)\ 



7r2 flTG {x) 

Number conservation / dxAh{x) — yields, 

64V25k 



Afig = 



(4.12) 



(4.13) 



Restoring the equations to dimensional form, the density 
profile of a strongly interacting gas becomes. 



n (x) 



(5>1 



n-TG {x) + An (x) , 



where 



n-TG (x) = y/Kn' 
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KX 
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Figure 5: (Color online) Density profiles of a ID trapped 
multi-component Fermi cloud at three interaction parameters 
5 = 0.1 (a), 5 = 1.0 (b), and 5 = 10 (c). The linear density 
and the coordinate are in units of the peak density n^F,K and 
Thomas- Fermi radius x%p^^ of an ideal gas, respectively. 



is the profile of a spinless Tonks-Girardeau gas, and the 
density variation, 
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KX 
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(4.16) 
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Figure 6: (Color online) Thomas-Fermi radius, in units of 
^TF.K, as a function of the interaction parameter, at several 
number of component as indicated. 




Figure 7: (Color online) Chemical potential as a function of 
the interaction parameter S. It is in units of -Ef.k = Nhw/ k,. 
Thin solid lines are the analytic results in the two limiting 
cases for k = 2, as described in Eqs. I|4.17p and l|4.2ip . 



Accordingly, the chemical potential takes the form, 



(k^ - 1) ft2 



Nhw 



64\/2kS'^ 1 
9^ '6 



(4.17) 

where the first term on the right hand side is again from 
the binding energy, while the second term corresponds 
to the chemical potential of a spinless Tonks-Girardeau 
gas of bound clusters. For later reference, we calculate 
the mean-square size of the cloud (x^) = / dxx'^n{x)/N 
using the strongly interacting density profile n{x)g^^, 
and find that. 



\^ /S»l - 2«;2°'»o+ 15^2^ 



Ar3/2 



(4.18) 



The density profile of a weakly interacting gas can be 
calculated in the same manner. The leading order con- 
tribution is simply the ideal gas result, riideaiix), and the 



11 



Hartree-Fock correction to the chemical potential gives 
rise to a density variation which is proportional to the 
interaction strength Qid- Explicitly, we find that, 



n (x) 



<5«1 



nideal (x) + An (x) , 



(4.19) 



where 



^ , ^ 2k(k-1) 
An (x) = — \ 



2/7r 
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The chemical potential is given by, 



^ 4^2^(^.-1) ^ 



The mean-square size of the cloud is found to be, 
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(4.21) 

(4.22) 



Fig. ^ plots numerical results for the density profiles 
at three interaction parameters and at different number 
of species as indicated. The linear density and the coordi- 
nate are in units of the peak density rt^^ ^ and Thomas- 
Fermi radius a;!^^^ of an ideal gas, respectively. With 
increasing interaction parameter S (from Figs. ([5^) to 
([5b)), the density profiles change from an ideal gas distri- 
bution to a strongly interacting Tonks-Girardeau profile. 
At the same interaction parameter, the density profiles 
become sharper and narrower as the number of species 
K increases. We display also the Thomas-Fermi radius 
of the cloud in Fig. ([6]) as a function of the interaction 
strength. It decreases monotonically from the ideal gas 
result Xj^p ^ to the Tonks-Girardeau prediction Xj^p ^/n 
as the interaction parameter increases, thus increasing 
the compressive force of the attractive interactions. 

Fig. ([7]) reports the dependence of the chemical po- 
tential on the interaction strength, at several numbers of 
species as indicated. Here, for clarity we have subtracted 
the binding energy part of the chemical potential. For 
K = 2, we show the asymptotic behavior given by Eqs. 
(|4T7|) and (|42T|) by thin solid lines. They fit very well 
with the numerical results, except at the intermediate 
interaction regime. 



B. Low-lying collective modes 

Experimentally, a useful way to characterize an inter- 
acting system is to measure its low-lying collective excita- 
tions of density oscillations. For a uniform gas, the low- 
lying collective excitations are simply the sound waves 
with energy 17(fc) = c |k| for a given momentum fc < kp^ 
which is gapless as A; ^ 0. In traps, however, the spec- 
trum becomes discrete, due to the finite cloud size of the 



gas that is of order x^^p ^. There is a minimum value of 
the momentum fcmin ^ I/^^tfk- Taking a sound velocity 
at the trap center c ~ {hn^p ^^f m){Tr / k) , one finds that. 



^(^min) 



1 



(4.23) 



comparable to the energy level of the harmonic trap. 

Since the charge degree of freedom of the gas falls into 
the Luttinger liquid universality class, quantitative cal- 
culations of the low-lying collective excitations in traps 
can be carried out based on the superfiuid hydrodynamic 
description of the dynamics of the ID Fermi gas [5ll.[60|. 
In such a description, the density n {x, t) and the velocity 
field v{x,t) satisfy the equation of continuity 



dn {x, t) d 
dt 9^ ' 

and the Euler equation 

dv d 



n {x, t) V (x, t)] = 0, 



dt dx 



Aihom ("-) + Vtrap {x) + -mv' 



(4.24) 



= 0. (4.25) 



We consider the fiuctuations of the density and the veloc- 
ity field about the equilibrium ground state , 5n (x, t) = 
n (x, t) — n(x) and 6v (x, t) = v (x, t) — v{x) = v (x, t), 
where n(x) and v{x) = are the equilibrium density 
profile and velocity field. Linearizing the hydrodynamic 
equations, one finds that (60| . 



hom(,n-) 



dn 



6n (x, t) 



^Sn(x t) = -—L— 
dt'^ ' m dx \ dx 

~ (4.26) 

The boundary condition requires that the current 
j{x,t) — n{x)Sv {x,t) should vanish identically at the 
Thomas-Fermi radius x = zLxpF- Considering the nth 
eigenmode with Sn (x, t) = Sn (x) exp [iujnt] and remov- 
ing the time-dependence, we end up with an eigenvalue 
problem, i.e., 



1 d 

m dx 



dx 



dn 



Sn (x) 



+ ujlSn{x) = 0. (4.27) 



We develop a powerful multi-series-expansion method 
to solve the above ID hydrodynamics equation, as out- 
lined in detail in the Appendix. The resulting low-lying 
collective mode can be classified by the number of nodes 
in its eigenfunction, i.e., the number index "n". The 
lowest two modes with n — 1,2 have very transparent 
physical meaning. These are respectively the dipole and 
breathing (compressional) modes, which can be excited 
separately by shifting the trap center or modulating the 
harmonic trapping frequency. The dipole mode is not af- 
fected by interactions according to Kohn's theorem, and 
has an invariant frequency precisely at wi = cj. 

The low-lying hydrodynamic modes of a two- 
component Fermi gas in the weak and strong coupling 
limits have been discussed analytically by Minguzzi [sij . 
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Figure 8: (Color online) Frequency corrections (with respect 
to the ideal gas value) of the lowest breathing modes Slob as a 
function of the dimensionless coupling parameter 5, at several 
number of species as indicated. 




Figure 9: (Color online) Interaction strength dependence of 
the frequency corrections Su)n = uj„ — nui of the low-lying 
collective modes of a ID three-component Fermi gas. Thin 
solid lines are the analytic sum-rule expansions for the weakly 
and strongly interacting limits, as described in Eqs. I|4.29p 
and l|4.30p . respectively. 



In both limits, the cloud behaves like a spinless ideal 
Fermi gas. Therefore, the frequency w„ of the mode "n" 
is fixed to nuj. This result applies to a multi-component 
ID Fermi gas as well. 

Fig. ^ shows the interaction dependence of frequen- 
cies of the breathing mode for a multi-component gas 
with different numbers of components. Here, to stress the 
role of interactions, the deviation of the mode frequency 
from its ideal result, Scub = ujb — has been plot- 
ted. As a function of the interaction parameter, a peak 
in Slub emerges at the intermediate interaction regime 
5 ~ 1. The peak value increases with increasing the 
number of components k. This peak is a clear signature 
of the cross-over from a regime of multi-particle clusters 
to a color quasi-superconductor. 



Fig. ([9]) shows the frequency correction Sujn = LOn — nuj 
of the low-lying modes of a three-component gas as a 
function of the interaction strength. The dipole mode 
frequency is always uj, as expected, while all other modes 
show a non-trivial peak structure, similar to that of the 
breathing mode. It is evident that the higher mode index 
"n" is, the larger frequency correction. 

As an alternative to the numerical solution of the ID 
hydrodynamic equation, the frequency of the breathing 
mode may be obtained by applying the compressibility 
sum-rule. As long as the breathing mode exhausts all 
the weights in the dynamic structure factor, the mode 
frequency can be calculated according to [isl. [60|. 
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(4.28) 



We have performed such calculations. The resulting fre- 
quency agrees extremely well with that from solving the 
ID hydrodynamic equation, with a relative difference less 
than 10""*. Such a good agreement gives strong support 
to the application of the sum-rule. To gain further insight 
of the breathing mode frequency, we use the sum-rule to 
obtain analytically the first order correction to the fre- 
quency for the weak and strong coupling limits. This 
can be done by substituting the expression for the mean- 
square size of the cloud in Eqs. I|4.18p and (|4.22p into 
the sum-rule formalism (|4.28p . Replacing aho = \/h/mLu 
and taking the derivative with respect to , we find that 



,5>1 



2w 1 



16\/2kS'« 1 
157r2 ~5 



for the strongly interacting limit, and 



,<5«1 



= 2uj 



1 ^ Tl. ^ 



(4.29) 



(4.30) 



for the weak coupling regime. Thin solid lines in Fig. 
((9| show the resulting analytic expansions for 5ujb for a 
three-component Fermi gas. These expansions describe 
very well the breathing mode frequency over a fairly large 
range of interaction strengths, but break down at the 
intermediate coupling regime 5^1. 



V. CONCLUSION 

In conclusion, we have investigated the properties of 
ID attractive multi-component Fermi gases, based on the 
Bethe ansatz exact solution, in a homogeneous environ- 
ment. This was extended to include a harmonic trap, by 
using the local density approximation. The equation of 
state of the system has been discussed in detail, as well as 
some dynamical quantities, including the sound velocity 
and low-lying collective modes. 

We have drawn attention to the formation of multi- 
body bound clusters, which are always present in a ID 
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attractive multi-component gas. These clusters are not 
well understood, as their description is beyond the mean- 
field theory. We have found that such multi-body clus- 
ters have significant impact on the equation of state and 
dynamic behavior of the systems. In particular, as the 
number of particles in the clusters increase, the gas turns 
out to be increasingly attractive compared to a sim- 
ple Cooper-pairing scenario. This is suggestive of the 
strongly interacting nature of the bound clusters. 

There is a close analogy between the proposed proper- 
ties of QCD at high density, and the results we calculate 
for a three-component Fermi gas in ID. In both cases 
there is a transition between discrete multi-particle clus- 
ters (nucleons) to a coherent quasi-superfiuid (colour su- 
perconductor), as the density increases. For the ID case, 
both the high and low density limits result in free particle 
dynamics, with different multiplicities. This gives simple 
collective mode behaviour in either limit. However, the 
transition region with 7 ~ 1 can show very strong evi- 
dence of a transition. This is due to an easily observed 
peak in the collective breathing mode frequency. 

Our results should be useful for the future experiments 
on ID multi-component atomic Fermi gases. An exam- 
ple of particular interest is three-component lithium gas, 
which has three broad Feshbach resonances. Our esti- 
mate of the relevant parameters suggests that an ultra- 
cold gas of ^Li atoms in an optical lattice, above a mag- 
netic field B = 1000 Gauss, can be nearly described by 
the current model. Our prediction of the density profiles 
and low-lying collective modes, provide a useful charac- 
terization of an interacting ID three-component Fermi 
gas. We expect this to be testable in a future experi- 
ment. 
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Appendix 

In this appendix, we outline the procedure of solving 
the ID hydrodynamic equation. 



m dx 



dx 



dn 



5n (x) 



+ ulSn{x) = 0, (5.1) 



with the boundary condition that the current j{x) oc 
dx{Sn (x) diihom[n'{x)]/dn{x)} must vanish at a; = zLxtf- 
It is convenient to introduce a function, 



^Mhom [n jx)] 
dn {x) 



Sn (x) 



(5.2) 



and rewrite the hydrodynamic equation in the form. 



d^/ dlii[n{x)] df fujm\'^ ^"^x 



X-rr; -H X- 

dx-^ 



dx dx 



(?-(x) 



/ = 0, (5.3) 



where c(x) = \n(x)d\i\^om\n{x)\l dn{x) jrn)^^'^ is the local 
sound velocity at position x. Let us introduce A (x) = 
ln[n(x)], A\ (x) = —lo^x/c^ {x), and Vm = LOm/^- Then 
the above equation becomes. 



dV dA df 



dx^ 



dx dx 



ulA, {x) / - 0. 



(5.4) 



Here the variable x £ [—xtf,+xtf]- Note that for our 
specific case, Ai {x) — dA{x)/dx, due to the local equi- 
librium condition. 

To solve the equation, one may wish to change x to a 
new variable y € [0, 1]. This can be done by two steps. 
First, we define / (x) = x^v (x), where Z = or 1 corre- 
sponds to the parity of the modes. Translating to v (x), 
we have (now x e [0, -|-xtf]), 



dx^ 



21 



dA\ dv 



dx j dx 



4i - I 



dA 

dx 



u = 0. (5.5) 



At the second step, we define y = (x/xtf) and change 
to the new variable y. After some straightforward calcu- 
lations, we obtain the following equation for v {y), 



2/(1 - 2/)^ + (1 



V) 



A 



V 



dA 
dy 



,dA 
dy 



^laM iy) 



dv 
dy 

, (5.6) 



where A = I + 1/2 and Ai{y) = Ai{y)/ {2y^^'^). In 
practise, the value of dA/dy can be calculated as follows. 



dA _ mLo^x'^p 



1 



nd^ihom[n]/dn 



(5.7) 

Note that in our case Ai (y) — dA/dy. Here, we multiply 
a factor of (1 — y)^ on both sides of Eq. (|5.6p . in order 
to remove the singularity of dA/dy and Ai {y) at point 

y = i. 

We develop a multi-series-expansion method to solve 
the eigenvalue problem Eq. I|5.6p . As the current van- 
ishes at the Thomas-Fermi boundary, the eigenfunction 
of Eq. (|5.6p should not be singular aX y — 1. As well, we 
require that the eigenfunction has to take a finite value 
at 2/ = 0. As we shall see, these two boundary conditions 
give rise to a set of the discrete spectrum, as well as a 
complete set of the eigenfunctions. 

To apply the boundary conditions, we divide the whole 
region [0, 1] into many pieces, for example, M parts, 
[0,1] = [2/0 = 0,2/1] U [2/1,2/2] U ... U [yM-i,yM = 1]. We 
look for the solution in the form, 



v{y) 



n=0 



(2/ - 2/i)" - if 2/ C [?/», 2/i+i] • (5.8) 
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Hence, the two boundary conditions translate to the re- 
quirement of a well-convergent series of {ai„} at both 
the starting region [yo = 0, yi] and the ending region 
[yM-i, yM = !]• The basic idea of solving the eigenvalue 
problem is then clear. We use the strategy of try and 
test. Given the parity I, we make an initial guess for 
i/m, and setup the series {ai„} at the starting region 
[yo = 0, yi], and then propagate it to the ending region of 
[um-i, yM = !]• If the series converges at y = 1, then we 
find a correct eigenvalue and eigenfunction of the prob- 
lem. Otherwise, we scan the value of i^m, until all the 
required low-lying eigenvalues are found. 

In greater detail, we apply the try and test strategy as 
follows. (A) At first, let us approximate, at each region 

[yi,yi+i], 



-^'~'^^y 
-(i-yfii {y) 



pq + piy + P2y'^ , (5.9) 



qa + qiy + hy'^ ■ (s.io) 



To make the expansion accurate, generally we take M ~ 
30. By introducing a new variable z = y — yi, at the 
region [yi,yi+i\ we can cast the Eq. I|5.6p into the form, 

E-.-^j + (g^^"') Ty + (g*^' ' " - 

(5.11) 

where the coefficients {ri \ , {pi\ and {qi\ can be calcu- 
lated directly from {pi} and {qi} in the program. We 
then substitute the solution v {z) = '}2^=o (^inz" into the 
above equation to obtain the iterative relation (without 



confusion, here we denote a„ = a.m for this region), 

_ {n + 1) (nri + po) 

an+2 — -7 —TT-, — T fln-l-l 

[n + 2){n+ 1) To 
[n {n- l)r2 + npi + go] 
(n + 2) (n + 1) ro 

[(n-l)(n-2)r3 + (n-l)p2+ Qi] 

(n + 2) {n+l)ro 
{{n-2)p3 + q2] 



-a-n-i 



{n + 2) (n+ l)ro' 



(5.12) 



We need to classify two cases, (i) In the starting region 
of yo = 0, we have tq = due to the boundary condition. 
Up to an overall irrelevant factor, we can set ag = 1 and 
then ai = —qo/po- (ii) In the other regions, aa and oi 
are determined by the continuous conditions as stated 
below. Once ao and ai are known, we could obtain all 
the values of a„ by Eq. (|5.12p since a_i — a_2 = 0. 
Usually it is already sufficiently accurate to keep n < 
'T-max = 16. (B) The series {ai„} at different regions 
are connected by the requirement that the function v (y) 
and its first derivation v' {y) should be continuous at the 
point {yi}, where the index i runs from 1 to M — 1. (C) 
In this way, we can finally obtain the series {ai„} at the 
region [?/m-1j 2/m]- We judge the convergence by checking 
whether the value of aA/,n„ax is sufficiently small or not. 

In practise, the above procedure of solving the ID hy- 
drodynamic equation of a multi-component Fermi gas is 
very efficient and accurate. It can be applied as well to 
other ID systems, such as the ID interacting Bose gas, 
and to the 3D systems in spherical harmonic traps. 
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